!interpolants
! E01AEF
! E01AEU
! E01AEV
! E01AEW
! E01AEX
! E01AEY
! E01AEZ

      SUBROUTINE E01AEF(M,XMIN,XMAX,X,Y,IP,N,ITMIN,ITMAX,A,WRK,LWRK,
     *                  IWRK,LIWRK,IFAIL)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C
C     *******************************************************
C
C     NPL ALGORITHMS LIBRARY ROUTINE PINTRP
C
C     CREATED 17 07 79.  UPDATED 14 05 80.  RELEASE 00/42
C
C     AUTHORS ... GERALD T. ANTHONY, MAURICE G. COX,
C                 J. GEOFFREY HAYES AND MICHAEL A. SINGER.
C     NATIONAL PHYSICAL LABORATORY, TEDDINGTON,
C     MIDDLESEX TW11 OLW, ENGLAND
C
C     *******************************************************
C
C     E01AEF. A ROUTINE, WITH CHECKS, WHICH DETERMINES AND
C     REFINES A POLYNOMIAL INTERPOLANT  Q(X)  TO DATA WHICH
C     MAY CONTAIN DERIVATIVES.
C
C     INPUT PARAMETERS
C        M        NUMBER OF DISTINCT X-VALUES
C        XMIN,
C        XMAX     LOWER AND UPPER ENDPOINTS OF INTERVAL
C        X        INDEPENDENT VARIABLE VALUES (DISTINCT)
C        Y        VALUES AND DERIVATIVES OF
C                    DEPENDENT VARIABLE.
C        IP       HIGHEST ORDER OF DERIVATIVE AT EACH X-VALUE.
C        N        NUMBER OF INTERPOLATING CONDITIONS.
C                    N = M + IP(1) + IP(2) + ... + IP(M).
C        ITMIN,
C        ITMAX    MINIMUM AND MAXIMUM NUMBER OF ITERATIONS TO BE
C                    PERFORMED.
C
C     OUTPUT PARAMETERS
C        A        CHEBYSHEV COEFFICIENTS OF  Q(X)
C
C     WORKSPACE (AND ASSOCIATED DIMENSION) PARAMETERS
C        WRK      REAL WORKSPACE ARRAY.  THE FIRST IMAX ELEMENTS
C                    CONTAIN, ON EXIT, PERFORMANCE INDICES FOR
C                    THE INTERPOLATING POLYNOMIAL, AND THE NEXT
C                    N  ELEMENTS THE COMPUTED RESIDUALS
C        LWRK     DIMENSION OF WRK. LWRK MUST BE AT LEAST
C                    7*N + 5*IMAX + M + 2, WHERE
C                    IMAX IS ONE MORE THAN THE LARGEST ELEMENT
C                    OF THE ARRAY IP.
C        IWRK     INTEGER WORKSPACE ARRAY.  ON EXIT,  IWRK(1)
C                    CONTAINS THE NUMBER OF ITERATIONS TAKEN
C        LIWRK    DIMENSION OF IWRK.  AT LEAST 2*M + 2.
C
C     FAILURE INDICATOR PARAMETER
C        IFAIL    FAILURE INDICATOR.
C                    0 - SUCCESSFUL TERMINATION.
C                    1 - AT LEAST ONE OF THE FOLLOWING CONDITIONS
C                        HAS BEEN VIOLATED -
C                           M AT LEAST 1,
C                           N = M + IP(1) + IP(2) + ... + IP(M),
C                           LWRK AT LEAST 7*N + 5*IMAX + M + 2,
C                           LIWRK AT LEAST 2*M + 2.
C                    2 - FOR SOME I, IP(I) IS LESS THAN 0.
C                    3 - AT LEAST ONE OF THE FOLLOWING CONDITIONS
C                        HAS BEEN VIOLATED -
C                           XMIN STRICTLY LESS THAN XMAX,
C                           FOR EACH I, X(I) MUST LIE IN THE
C                              INTERVAL XMIN TO XMAX,
C                           THE X-VALUES MUST ALL BE DISTINCT
C                    4 - NOT ALL PERFORMANCE INDICES LESS THAN
C                        ONE, BUT ITMAX ITERATIONS PERFORMED,
C                    5 - COMPUTATION TERMINATED BECAUSE
C                        ITERATIONS DIVERGING.
C
C
C     CHECK AND SET ITERATION LIMITS
C
C     .. Parameters ..
      CHARACTER*6       SRNAME
      PARAMETER         (SRNAME='E01AEF')
C     .. Scalar Arguments ..
      DOUBLE PRECISION  XMAX, XMIN
      INTEGER           IFAIL, ITMAX, ITMIN, LIWRK, LWRK, M, N
C     .. Array Arguments ..
      DOUBLE PRECISION  A(N), WRK(LWRK), X(M), Y(N)
      INTEGER           IP(M), IWRK(LIWRK)
C     .. Local Scalars ..
      DOUBLE PRECISION  XI
      INTEGER           I, IERROR, IMAX, IP1, ITMAX1, ITMIN1, J, NN
C     .. Local Arrays ..
      CHARACTER*1       P01REC(1)
C     .. External Functions ..
      INTEGER           P01ABF
      EXTERNAL          P01ABF
C     .. External Subroutines ..
      EXTERNAL          E01AEW
C     .. Executable Statements ..
!       write(*,*) 'E01AEF 0 N,M=',N,M
      ITMIN1 = ITMIN
      ITMAX1 = ITMAX
      IF (ITMIN1.LE.0) ITMIN1 = 2
      IF (ITMAX1.LE.0) ITMAX1 = 10
      IF (ITMAX1.LT.ITMIN1) ITMAX1 = ITMIN1
C
C     FIRST SET OF CHECKS ON INPUT DATA
C
      IERROR = 1
C
      IF (M.LT.1) GO TO 100
C
C     DETERMINE  IMAX = MAX(IP(I) + 1),  AND  NN = SUM
C     OF VALUES OF  (IP(I) + 1)  FOR COMPARISON WITH  N
C
      IMAX = 0
      NN = M
      DO 20 I = 1, M
         IF (IP(I).GT.IMAX) IMAX = IP(I)
         NN = NN + IP(I)
   20 CONTINUE
      IMAX = IMAX + 1
C
C     REMAINDER OF FIRST SET OF INPUT DATA CHECKS
C
      IF (NN.NE.N) GO TO 100
      IF (LWRK.LT.7*N+5*IMAX+M+2) GO TO 100
      IF (LIWRK.LT.2*M+2) GO TO 100
C
C     SECOND SET OF CHECKS ON INPUT DATA
C
      IERROR = 2
      DO 40 I = 1, M
         IF (IP(I).LT.0) GO TO 100
   40 CONTINUE
C
C     THIRD SET OF CHECKS ON INPUT DATA
C
      IERROR = 3
!       IF (XMIN.GE.XMAX) write(*,*) 'E01AEF 3 1'
      IF (XMIN.GE.XMAX) GO TO 100
      DO 80 I = 1, M
         XI = X(I)
!          IF (XI.LT.XMIN .OR. XI.GT.XMAX) write(*,*) 'E01AEF 3 2'
         IF (XI.LT.XMIN .OR. XI.GT.XMAX) GO TO 100
         IF (I.EQ.M) GO TO 80
         IP1 = I + 1
         DO 60 J = IP1, M
!             IF (XI.EQ.X(J)) write(*,*) 'E01AEF 3 3 I, J, XI, X(J)', I, J, XI, X(J)
            IF (XI.EQ.X(J)) GO TO 100
   60    CONTINUE
   80 CONTINUE
C
C     END OF CHECKS ON INPUT DATA
C
      IERROR = 0
C
C     INDICATE THAT ZEROIZING POLYNOMIAL  Q0(X)  IS
C     NOT REQUIRED
C
      IWRK(1) = 0
      CALL E01AEW(M,XMIN,XMAX,X,Y,IP,N,N+1,ITMIN1,ITMAX1,A,WRK,WRK,LWRK,
     *            IWRK,LIWRK,IERROR)
C
C     SET APPROPRIATE VALUE FOR FAILURE INDICATOR
C
      IF (IERROR.NE.0) IERROR = IERROR + 3
  100 IFAIL = P01ABF(IFAIL,IERROR,SRNAME,0,P01REC)
      RETURN
C
C     END E01AEF
C
      END

      SUBROUTINE E01AEU(M,X,IP,NP1,B,W)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C     MARK 13 REVISED. USE OF MARK 12 X02 FUNCTIONS (APR 1988).
C
C     *******************************************************
C
C     NPL ALGORITHMS LIBRARY ROUTINE Q0POLY
C
C     CREATED 02 05 80.  UPDATED 13 05 80.  RELEASE 00/08
C
C     AUTHORS ... GERALD T. ANTHONY, MAURICE G. COX
C                 J. GEOFFREY HAYES AND MICHAEL A. SINGER.
C     NATIONAL PHYSICAL LABORATORY, TEDDINGTON,
C     MIDDLESEX TW11 OLW, ENGLAND.
C
C     *******************************************************
C
C     E01AEU.  AN ALGORITHM TO DETERMINE THE CHEBYSHEV SERIES
C     REPRESENTATION OF A ZEROIZING POLYNOMIAL  Q0(X),
C     I.E. A POLYNOMIAL WHICH TAKES ON ZERO VALUES (AND
C     POSSIBLY ZERO DERIVATIVE VALUES) AT SPECIFIED POINTS
C
C     INPUT PARAMETERS
C        M        NUMBER OF DISTINCT X-VALUES.
C        X        INDEPENDENT VARIABLE VALUES,
C                    NORMALIZED TO  (-1, 1)
C        IP       HIGHEST ORDER OF DERIVATIVE AT EACH X-VALUE
C        NP1      N + 1,  WHERE  N = NUMBER OF ZEROS (INCLUDING
C                    THOSE OF DERIVATIVES) TO BE TAKEN ON BY
C                    Q0(X).  N = M + IP(1) + IP(2) + ... + IP(M).
C
C     OUTPUT PARAMETERS
C        B        CHEBYSHEV COEFFICIENTS OF  Q0(X)
C
C     WORKSPACE PARAMETERS
C        W        WORKSPACE
C
C     .. Scalar Arguments ..
      INTEGER           M, NP1
C     .. Array Arguments ..
      DOUBLE PRECISION  B(NP1), W(NP1), X(M)
      INTEGER           IP(M)
C     .. Local Scalars ..
      DOUBLE PRECISION  AI, ANBIG, CT, EPS, FACTOR, ONE, OVFL, PI, RI,
     *                  SFAC, SXTEEN, SXTNTH, TEST, TWO, UNFL, XI,
     *                  XTREMA, ZERO
      INTEGER           I, I2, IFAIL, IP1, K, L, N, NU
C     .. External Functions ..
      DOUBLE PRECISION  X01AAF, X02AJF, X02AMF
      EXTERNAL          X01AAF, X02AJF, X02AMF
C     .. External Subroutines ..
      EXTERNAL          E02AFF
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, LOG, SIN
C     .. Data statements ..
      DATA              ZERO, SXTNTH, ONE, TWO, SXTEEN/0.0D+0,
     *                  0.0625D+0, 1.0D+0, 2.0D+0, 16.0D+0/
C     .. Executable Statements ..
C
C     EVALUATE  Q0(X)  AT THE EXTREMA OF THE CHEBYSHEV POLYNOMIAL
C     OF DEGREE  N,  EMPLOYING DYNAMIC SCALING TO AVOID THE
C     POSSIBILITY OF OVERFLOW OR UNDERFLOW
C
      OVFL = SXTNTH/X02AMF()
      UNFL = SXTNTH*OVFL
      EPS = X02AJF()
      PI = X01AAF(ZERO)
      N = NP1 - 1
      FACTOR = 2*N
      FACTOR = PI/FACTOR
      DO 20 I = 1, NP1
         W(I) = ONE
   20 CONTINUE
      DO 120 K = 1, M
         IP1 = IP(K) + 1
         XI = X(K)
         DO 100 L = 1, IP1
            ANBIG = ZERO
            I2 = N + 2
            DO 40 I = 1, NP1
               I2 = I2 - 2
               RI = I2
               XTREMA = SIN(FACTOR*RI)
               AI = W(I)*(XTREMA-XI)
               W(I) = AI
               IF (ABS(AI).GT.ANBIG) ANBIG = ABS(AI)
   40       CONTINUE
   60       SFAC = ONE
            IF (ANBIG.GT.OVFL) SFAC = SXTNTH
            IF (ANBIG.LT.UNFL) SFAC = SXTEEN
            IF (SFAC.EQ.ONE) GO TO 100
            ANBIG = ANBIG*SFAC
            DO 80 I = 1, NP1
               W(I) = W(I)*SFAC
   80       CONTINUE
            GO TO 60
  100    CONTINUE
  120 CONTINUE
      CT = OVFL
      DO 140 I = 1, NP1
         TEST = ABS(W(I))/ANBIG
         IF (TEST.LE.EPS) W(I) = ZERO
         IF (TEST.GT.EPS .AND. TEST.LT.CT) CT = TEST
  140 CONTINUE
      CT = CT*ANBIG
      SFAC = ONE
  160 IF (CT.LT.ONE) GO TO 180
      CT = CT*SXTNTH
      SFAC = SFAC*SXTNTH
      GO TO 160
  180 DO 200 I = 1, NP1
         W(I) = W(I)*SFAC
  200 CONTINUE
C
C     DETERMINE THE CHEBYSHEV REPRESENTATION OF  Q0(X)
C
      CALL E02AFF(NP1,W,B,IFAIL)
C
C     SET THE LEADING COEFFICIENT OF  Q0(X)  TO AN
C     EXACT POWER OF  2
C
      AI = B(NP1)
      CT = LOG(ABS(AI))/LOG(TWO)
      NU = CT
      SFAC = TWO**NU
      B(NP1) = SFAC
      SFAC = SFAC/AI
      DO 220 I = 1, N
         B(I) = B(I)*SFAC
  220 CONTINUE
C
      RETURN
C
C     END E01AEU
C
      END


      SUBROUTINE E01AEV(WITHQ0,WITHPI,M,XMIN,XMAX,X,N,Y,IP,IMAX,A,LA,IT,
     *                  RNMBST,RNM,IMPROV,ADIF,RES,PMAX,PINDEX)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C     MARK 13 REVISED. USE OF MARK 12 X02 FUNCTIONS (APR 1988).
C
C     *******************************************************
C
C     NPL ALGORITHMS LIBRARY ROUTINE PRESID
C
C     CREATED 18 02 80.  UPDATED 14 05 80.  RELEASE 00/28
C
C     AUTHOR ... MAURICE G. COX.
C     NATIONAL PHYSICAL LABORATORY, TEDDINGTON,
C     MIDDLESEX TW11 OLW, ENGLAND
C
C     *******************************************************
C
C     E01AEV.  FORMS PERFORMANCE INDICES AND RESIDUALS
C     FOR A POLYNOMIAL APPROXIMATION  P(X)  TO A SET OF DATA
C     WHICH MAY CONTAIN DERIVATIVE VALUES.  ALSO INDICATES
C     WHETHER THE POLYNOMIAL DEFINED BY THE SPECIFIED
C     COEFFICIENTS IS A BETTER APPROXIMATION THAN THE BEST
C     SO FAR.
C
C     INPUT PARAMETERS
C        WITHQ0   TRUE IF ZEROIZING POLYNOMIAL, ELSE FALSE
C        WITHPI   IF TRUE, PERFORMANCE INDICES PRODUCED UNDER
C                    ALL CIRCUMSTANCES.  OTHERWISE, PRODUCED
C                    ONLY IF  P(X)  IS AN IMPROVEMENT
C        M        NUMBER OF X-VALUES.  ALL DISTINCT
C        XMIN,
C        XMAX     LOWER AND UPPER ENDPOINTS OF INTERVAL
C        X        X-VALUES.  NORMALIZED TO  (-1, 1)
C        N        NUMBER OF Y-VALUES
C        Y        VALUES AND DERIVATIVES OF DEPENDENT VARIABLE
C        IP       HIGHEST ORDER OF DERIVATIVE AT EACH X-VALUE
C        IMAX     ONE GREATER THAN LARGEST VALUE OF  IP
C        A        CHEBYSHEV COEFFICIENTS OF  P(X)
C        LA       DIMENSION OF  A  AND  ADIF.
C                    .GE. N + 1 IF WITHQ0 IS TRUE,
C                    .GE. N     OTHERWISE
C        IT       ITERATION NUMBER
C
C     INPUT/OUTPUT PARAMETERS
C        RNMBST   2-NORMS OF RESIDUALS CORRESPONDING TO THE
C                    BEST POLYNOMIAL SO FAR AND ITS DERIVATIVES
C
C     OUTPUT PARAMETERS
C        RNM      2-NORMS OF RESIDUALS CORRESPONDING TO
C                    P(X)  AND ITS DERIVATIVES
C        IMPROV   TRUE IF  P(X)  IS AN IMPROVEMENT, ELSE FALSE
C        ADIF     CHEBYSHEV COEFFICIENTS OF  (IMAX - 1)-ST
C                    DERIVATIVE OF  P(X)
C        RES      RESIDUALS CORRESPONDING TO Y-VALUES
C      * PMAX     LARGEST PERFORMANCE INDEX
C      * PINDEX   PERFORMANCE INDICES
C
C     NOTE.  THE PARAMETERS MARKED  *  ARE PROVIDED ONLY
C            IF  IMPROV  IS TRUE
C
C     KEY LOCAL VARIABLES
C        ASUMAX   FOR CURRENT VALUE OF  L,  MAXIMUM VALUE
C                    OVER  J = 0, 1, ..., L  OF SUM OF MODULI
C                    OF CHEBYSHEV COEFFICIENTS OF DERIVATIVE
C                    OF ORDER  J  OF  P(X)
C        EPS      RELATIVE MACHINE PRECISION
C        IY       LOCATION IN  Y  OF CURRENT SPECIFIED
C                    DERIVATIVE VALUE OF ORDER  L
C        NL       NUMBER OF SPECIFIED DERIVATIVE VALUES
C                    OF ORDER  L
C        NTERMS   NUMBER OF TERMS IN CHEBYSHEV SERIES
C                    REPRESENTATION OF CURRENT DERIVATIVE
C                    OF  P(X)
C        RES2NM   2-NORM OF THE RESIDUALS CORRESPONDING TO THE
C                    SPECIFIED DERIVATIVE VALUES OF ORDER  L - 1
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION  PMAX, XMAX, XMIN
      INTEGER           IMAX, IT, LA, M, N
      LOGICAL           IMPROV, WITHPI, WITHQ0
C     .. Array Arguments ..
      DOUBLE PRECISION  A(LA), ADIF(LA), PINDEX(IMAX), RES(N),
     *                  RNM(IMAX), RNMBST(IMAX), X(M), Y(N)
      INTEGER           IP(M)
C     .. Local Scalars ..
      DOUBLE PRECISION  ABSRES, ASUM, ASUMAX, EPS, HALF, MLTPLR, ONE, P,
     *                  PLARGE, RESMAX, RSCALE, T, ZERO
      INTEGER           I, IA, IY, L, LM1, NL, NTERMS
      LOGICAL           IMP
C     .. External Functions ..
      DOUBLE PRECISION  X02AJF
      EXTERNAL          X02AJF
C     .. External Subroutines ..
      EXTERNAL          E02AHZ, E02AKZ
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, SQRT
C     .. Data statements ..
      DATA              HALF, ZERO, ONE, MLTPLR/0.5D+0, 0.0D+0, 1.0D+0,
     *                  8.0D+0/
C     .. Executable Statements ..
      PMAX = ZERO
      EPS = X02AJF()
      NTERMS = N
      IF (WITHQ0) NTERMS = NTERMS + 1
      ASUMAX = ZERO
      DO 20 I = 1, NTERMS
         ADIF(I) = A(I)
   20 CONTINUE
      NTERMS = NTERMS + 1
      DO 120 L = 1, IMAX
         NTERMS = NTERMS - 1
C
C        DETERMINE SUM OF MODULI  ASUM  OF CHEBYSHEV COEFFICIENTS
C        OF DERIVATIVE OF ORDER  L - 1  OF  P(X),  AND UPDATE
C        ASUMAX
C
         ASUM = HALF*ABS(ADIF(1))
         IF (NTERMS.EQ.1) GO TO 60
         DO 40 IA = 2, NTERMS
            ASUM = ASUM + ABS(ADIF(IA))
   40    CONTINUE
   60    IF (ASUM.GT.ASUMAX) ASUMAX = ASUM
C
C        PINDEX(L)  IS USED TEMPORARILY TO HOLD THE
C        VALUE OF  ASUMAX
C
         PINDEX(L) = ASUMAX
         IY = L
         NL = 0
C
C        TO REDUCE THE POSSIBILITY OF UNDERFLOW AND OVERFLOW,
C        COMPUTE  RES2NM  AS  RESMAX*SQRT(RSCALE),  WHERE
C        RESMAX  AND  RSCALE  ARE UPDATED FOR EACH RESIDUAL
C        CORRESPONDING TO A SPECIFIED DERIVATIVE VALUE OF
C        ORDER  L - 1.  AT ANY STAGE,  RESMAX  HOLDS THE
C        MODULUS OF THE RESIDUAL OF MAXIMUM MAGNITUDE SO FAR,
C        AND  RSCALE  THE SUM SO FAR OF THE SQUARES OF THE
C        RESIDUALS, EACH SCALED BY  RESMAX.
C
         RESMAX = ZERO
         RSCALE = ONE
         DO 100 I = 1, M
C
C           SKIP IF NO DERIVATIVE VALUE OF ORDER  L - 1  IS
C           SPECIFIED AT  I-TH  X-VALUE
C
            IF (IP(I)+1.LT.L) GO TO 80
            NL = NL + 1
C
C           EVALUATE  P,  THE  (L - 1)-ST  DERIVATIVE OF  P(X)
C           AT  X = X(I)
C
            CALL E02AKZ(NTERMS,ADIF,1,LA,X(I),P)
C
C           SAVE RESIDUAL CORRESPONDING TO THIS VALUE
C
            IF (WITHQ0) RES(IY) = -P
            IF ( .NOT. WITHQ0) RES(IY) = Y(IY) - P
C
C           UPDATE  RESMAX  AND  RSCALE
C
            ABSRES = ABS(RES(IY))
            IF (ABSRES.NE.ZERO .AND. ABSRES.LE.RESMAX) RSCALE = RSCALE +
     *          (ABSRES/RESMAX)**2
            IF (ABSRES.NE.ZERO .AND. ABSRES.GT.RESMAX)
     *          RSCALE = RSCALE*(RESMAX/ABSRES)**2 + ONE
            IF (ABSRES.NE.ZERO .AND. ABSRES.GT.RESMAX) RESMAX = ABSRES
   80       IY = IY + IP(I) + 1
  100    CONTINUE
         RNM(L) = RESMAX*SQRT(RSCALE)
C
C        FORM COEFFICIENTS IN CHEBYSHEV SERIES REPRESENTATION
C        OF  L-TH  DERIVATIVE OF  P(X)  FROM THOSE OF
C        (L - 1)-ST DERIVATIVE
C
         IF (L.LT.IMAX) CALL E02AHZ(NTERMS,XMIN,XMAX,ADIF,1,LA,T,ADIF,1,
     *                              LA)
  120 CONTINUE
C
C     IF NOT ON ZERO-TH ITERATION,
C     DETECT WHETHER THERE HAS BEEN AN IMPROVEMENT ...
C
      IMP = (IT.EQ.0)
      IF (IMP) GO TO 160
      DO 140 L = 1, IMAX
         IF (RNM(L).LT.RNMBST(L)) IMP = .TRUE.
  140 CONTINUE
      IF ( .NOT. (IMP .OR. WITHPI)) GO TO 220
C
C     ... AND, IF SO, OR IF ZERO-TH ITERATION, OR IF
C     THEY ARE REQUIRED ANYWAY, FORM THE PERFORMANCE
C     INDICES CORRESPONDING TO THE IMPROVED
C     APPROXIMATION
C
  160 PLARGE = ZERO
      DO 200 L = 1, IMAX
C
C        NL  IS THE NUMBER OF DERIVATIVE VALUES OF ORDER  L - 1
C
         NL = 0
         LM1 = L - 1
         DO 180 I = 1, M
            IF (IP(I).GE.LM1) NL = NL + 1
  180    CONTINUE
         T = NL
C
C        RETRIEVE THE VALUE OF ASUMAX
C
         ASUMAX = PINDEX(L)
         IF (ASUMAX.NE.ZERO) PINDEX(L) = RNM(L)
     *       /(MLTPLR*EPS*ASUMAX*SQRT(T))
         IF (ASUMAX.EQ.ZERO) PINDEX(L) = ZERO
         IF (PINDEX(L).GT.PLARGE) PLARGE = PINDEX(L)
  200 CONTINUE
      PMAX = PLARGE
  220 IMPROV = IMP
      RETURN
C
C     END E01AEV
C
      END


      SUBROUTINE E01AEW(M,XMIN,XMAX,X,Y,IP,N,NP1,ITMIN,ITMAX,A,B,WRK,
     *                  LWRK,IWRK,LIWRK,IFAIL)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C
C     *******************************************************
C
C     NPL ALGORITHMS LIBRARY ROUTINE PNTRPA
C
C     CREATED 20 12 79.  UPDATED 13 05 80.  RELEASE 00/47
C
C     AUTHORS ... GERALD T. ANTHONY, MAURICE G. COX,
C                 J. GEOFFREY HAYES AND MICHAEL A. SINGER.
C     NATIONAL PHYSICAL LABORATORY, TEDDINGTON,
C     MIDDLESEX TW11 OLW, ENGLAND
C
C     *******************************************************
C
C     E01AEW. A ROUTINE, WITHOUT CHECKS, WHICH DETERMINES AND
C     REFINES A POLYNOMIAL INTERPOLANT  Q(X)  TO DATA WHICH
C     MAY CONTAIN DERIVATIVES, AND AN ASSOCIATED ZEROIZING
C     POLYNOMIAL  Q0(X).
C
C        M        NUMBER OF DISTINCT X-VALUES
C        XMIN,
C        XMAX     LOWER AND UPPER ENDPOINTS OF INTERVAL
C        X        INDEPENDENT VARIABLE VALUES (DISTINCT)
C        Y        VALUES AND DERIVATIVES OF
C                    DEPENDENT VARIABLE.
C        IP       HIGHEST ORDER OF DERIVATIVE AT EACH X-VALUE.
C        N        NUMBER OF INTERPOLATING CONDITIONS.
C                    N = M + IP(1) + IP(2) + ... + IP(M).
C        NP1      VALUE OF  N + 1
C        ITMIN,
C        ITMAX    MINIMUM AND MAXIMUM NUMBER OF ITERATIONS TO BE
C                    PERFORMED.
C        IWRK(1)  SEE WORKSPACE (AND ASSOCIATED
C                    DIMENSION) PARAMETERS
C
C     OUTPUT PARAMETERS
C        A        CHEBYSHEV COEFFICIENTS OF  Q(X)
C        B        CHEBYSHEV COEFFICIENTS OF  Q0(X)
C
C     WORKSPACE (AND ASSOCIATED DIMENSION) PARAMETERS
C        WRK      REAL WORKSPACE ARRAY.  THE FIRST IMAX ELEMENTS
C                    CONTAIN, ON EXIT, PERFORMANCE INDICES FOR
C                    THE INTERPOLATING POLYNOMIAL, AND THE NEXT
C                    N  ELEMENTS THE COMPUTED RESIDUALS
C        LWRK     DIMENSION OF WRK. LWRK MUST BE AT LEAST
C                    8*N + 5*IMAX + M + 5, WHERE
C                    IMAX IS ONE MORE THAN THE LARGEST ELEMENT
C                    OF THE ARRAY IP.
C        IWRK     INTEGER WORKSPACE ARRAY.
C                    THE FIRST ELEMENT OF THIS ARRAY IS USED
C                    AS AN INPUT PARAMETER (WHICH IS DESTROYED
C                    ON EXIT).  THE ZEROIZING POLYNOMIAL  Q0(X)
C                    IS CONSTRUCTED OR NOT ACCORDING TO WHETHER
C                    IWRK(1)  IS NON-ZERO OR ZERO.
C        LIWRK    DIMENSION OF IWRK.  AT LEAST 2*M + 2.
C
C     FAILURE INDICATOR PARAMETER
C        IFAIL    FAILURE INDICATOR
C                    0 - SUCCESSFUL TERMINATION
C                    1 - ITERATION LIMIT IN DERIVING  Q(X)
C                    2 - DIVERGENT ITERATION IN DERIVING  Q(X)
C                    3 - ITERATION LIMIT IN DERIVING  Q0(X)
C                    4 - DIVERGENT ITERATION IN DERIVING  Q0(X)
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION  XMAX, XMIN
      INTEGER           IFAIL, ITMAX, ITMIN, LIWRK, LWRK, M, N, NP1
C     .. Array Arguments ..
      DOUBLE PRECISION  A(N), B(NP1), WRK(LWRK), X(M), Y(N)
      INTEGER           IP(M), IWRK(LIWRK)
C     .. Local Scalars ..
      DOUBLE PRECISION  ONE, PMAX, TWO, ZERO
      INTEGER           I, IADIF, IATRL, IBDIF, IBTRL, IC, ID, IDA, IDB,
     *                  IERROR, IFTAU, ILOCX, ILOCY, IMAX, INITQ,
     *                  INITQ0, IPIQ, IPIQ0, IPTRL, IRES, IRNM, IRTRNM,
     *                  IW, IX
      LOGICAL           WITHQ0
C     .. External Subroutines ..
      EXTERNAL          E01AEY, E02AKY
C     .. Data statements ..
      DATA              ZERO, ONE, TWO/0.0D+0, 1.0D+0, 2.0D+0/
C     .. Executable Statements ..
      IERROR = 0
      IMAX = 0
      DO 20 I = 1, M
         IF (IP(I).GT.IMAX) IMAX = IP(I)
   20 CONTINUE
      IMAX = IMAX + 1
C
C     INDICATE WHETHER  Q0(X)  IS REQUIRED
C
      WITHQ0 = (IWRK(1).NE.0)
C
C     TREAT THE CASE N = 1 SEPARATELY.
C
      IF (N.NE.1) GO TO 60
      A(1) = TWO*Y(1)
      IF ( .NOT. WITHQ0) GO TO 40
      CALL E02AKY(XMIN,XMAX,X(1),B(1))
      B(1) = -TWO*B(1)
      B(2) = ONE
C
C     SET TO ZERO THE NUMBERS OF ITERATIONS TAKEN, THE (SOLE)
C     RESIDUAL AND THE VALUES OF THE PERFORMANCE INDICES
C     IN THIS SPECIAL CASE, AND FINISH
C
   40 IWRK(1) = 0
      IF (WITHQ0) IWRK(2) = 0
      WRK(1) = ZERO
      WRK(2) = ZERO
      IF (WITHQ0) WRK(3) = ZERO
C
C     TRANSFORM THE  X*S  TO  (-1, 1)
C
   60 IX = 2*IMAX + N
      DO 80 I = 1, M
         IX = IX + 1
         CALL E02AKY(XMIN,XMAX,X(I),WRK(IX))
   80 CONTINUE
      IF (N.EQ.1) GO TO 120
      IF ( .NOT. WITHQ0) GO TO 100
C
C     WORKSPACE ALLOCATION FOR CALL TO  E01AEY  ...
C
      IRES = IMAX + 1
      IPIQ0 = IRES + N
      IX = IPIQ0 + IMAX
      IBTRL = IX + M
      IPTRL = IBTRL + NP1
      IFTAU = IPTRL + IMAX
      IC = IFTAU + N
      ID = IC + N
      IW = ID + NP1
      IBDIF = IW + NP1
      IDB = IBDIF + NP1
      IRNM = IDB + NP1
      IRTRNM = IRNM + IMAX
      INITQ0 = 2
      ILOCX = INITQ0 + 1
      ILOCY = ILOCX + M
C
C     ... TO DETERMINE  Q0(X)
C
      CALL E01AEY(.TRUE.,M,WRK(IX),XMIN,XMAX,Y,IP,IMAX,N,NP1,ITMIN,
     *            ITMAX,B,NP1,WRK(IRES),PMAX,WRK(IPIQ0),IWRK(INITQ0)
     *            ,WRK(IBTRL),WRK(IPTRL),WRK(IFTAU),WRK(IC),WRK(ID)
     *            ,WRK(IW),WRK(IBDIF),WRK(IDB),WRK(IRNM),WRK(IRTRNM)
     *            ,IWRK(ILOCX),IWRK(ILOCY),IERROR)
C
C     RE-CODE FAILURE INDICATOR
C
      IF (IERROR.NE.0) IERROR = IERROR + 2
      IF (IERROR.NE.0) GO TO 120
C
C     WORKSPACE ALLOCATION FOR CALL TO  E01AEY  ...
C
  100 IPIQ = 1
      IRES = IPIQ + IMAX
      IX = IRES + N + IMAX
      IATRL = IX + M
      IPTRL = IATRL + N
      IFTAU = IPTRL + IMAX
      IC = IFTAU + N
      ID = IC + N
      IW = ID + N
      IADIF = IW
      IDA = IADIF + NP1
      IRNM = IDA + NP1
      IRTRNM = IRNM + IMAX
      INITQ = 1
      ILOCX = INITQ + 2
      ILOCY = ILOCX + M
C
C     ... TO DETERMINE  Q(X)
C
      CALL E01AEY(.FALSE.,M,WRK(IX),XMIN,XMAX,Y,IP,IMAX,N,NP1,ITMIN,
     *            ITMAX,A,N,WRK(IRES),PMAX,WRK(IPIQ),IWRK(INITQ)
     *            ,WRK(IATRL),WRK(IPTRL),WRK(IFTAU),WRK(IC),WRK(ID)
     *            ,WRK(IW),WRK(IADIF),WRK(IDA),WRK(IRNM),WRK(IRTRNM)
     *            ,IWRK(ILOCX),IWRK(ILOCY),IERROR)
C
  120 IFAIL = IERROR
      RETURN
C
C     END E01AEW
C
      END


      SUBROUTINE E01AEX(M,X,IP,N,LOCX,C,NC,XNEW,IXNEXT,YNEW,NORDP1,CNEW,
     *                  D)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C
C     *******************************************************
C
C     NPL ALGORITHMS LIBRARY ROUTINE DIVDIF
C
C     CREATED 17 07 79.  UPDATED 14 05 80.  RELEASE 00/08
C
C     AUTHORS ... MAURICE G. COX AND MICHAEL A. SINGER.
C     NATIONAL PHYSICAL LABORATORY, TEDDINGTON,
C     MIDDLESEX TW11 OLW, ENGLAND
C
C     *******************************************************
C
C     E01AEX.  AN ALGORITHM TO DETERMINE THE NEXT COEFFICIENT
C     IN THE NEWTON FORM OF AN INTERPOLATING POLYNOMIAL
C
C     INPUT PARAMETERS
C        M        NUMBER OF DISTINCT X-VALUES.
C        X        INDEPENDENT VARIABLE VALUES,
C                    NORMALIZED TO  (-1, 1)
C        IP       HIGHEST ORDER OF DERIVATIVE AT EACH X-VALUE
C        N        NUMBER OF INTERPOLATING CONDITIONS.
C                    N = M + IP(1) + IP(2) + ... + IP(M).
C        LOCX     POINTERS TO X-VALUES IN CONSTRUCTING
C                    NEWTON FORM OF POLYNOMIAL
C        C        NEWTON COEFFICIENTS DETERMINED SO FAR
C        NC       NUMBER OF NEWTON COEFFICIENTS DETERMINED SO FAR
C        XNEW     ELEMENT OF  X  ASSOCIATED WITH NEW
C                    NEWTON COEFFICIENT
C        IXNEXT   NUMBER OF X-VALUES SO FAR INCORPORATED
C                    (INCLUDING  XNEW)
C        YNEW     SCALED DERIVATIVE VALUE CORRESPONDING TO
C                    XNEW  AND  NORDP1
C        NORDP1   ONE PLUS ORDER OF DERIVATIVE
C                    ASSOCIATED WITH  YNEW
C
C     INPUT/OUTPUT PARAMETERS
C        D        ELEMENTS IN PREVIOUS, AND THEN NEW, UPWARD
C                    SLOPING DIAGONAL OF DIVIDED DIFFERENCE TABLE
C
C     OUTPUT PARAMETERS
C        CNEW     NEW NEWTON COEFFICIENT GENERATED
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION  CNEW, XNEW, YNEW
      INTEGER           IXNEXT, M, N, NC, NORDP1
C     .. Array Arguments ..
      DOUBLE PRECISION  C(N), D(N), X(M)
      INTEGER           IP(M), LOCX(M)
C     .. Local Scalars ..
      DOUBLE PRECISION  DIF
      INTEGER           IC, IS, IX, K, LOCXI
C     .. Executable Statements ..
      IC = NC - NORDP1 + 1
      D(1) = YNEW
      IF (IXNEXT.EQ.1) GO TO 60
      IS = 0
      IX = 0
      DO 40 K = 1, IC
         IF (K.LE.IS) GO TO 20
         IX = IX + 1
         LOCXI = LOCX(IX)
         IS = IS + IP(LOCXI) + 1
         DIF = X(LOCXI) - XNEW
   20    IF (NORDP1.EQ.1) D(K+1) = (C(K)-D(K))/DIF
         IF (NORDP1.GT.1) D(K+1) = (D(K+1)-D(K))/DIF
   40 CONTINUE
   60 CNEW = D(IC+1)
      RETURN
C
C     END E01AEX
C
      END

      SUBROUTINE E01AEY(WITHQ0,M,X,XMIN,XMAX,Y,IP,IMAX,N,NP1,ITMIN,
     *                  ITMAX,A,LA,RES,PMAX,PINDEX,NIT,ATRIAL,PTRIAL,
     *                  FTAU,C,D,W,ADIF,DA,RNM,RTRLNM,LOCX,LOCY,IFAIL)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C
C     *******************************************************
C
C     NPL ALGORITHMS LIBRARY ROUTINE REFH
C
C     CREATED 17 07 79.  UPDATED 14 05 80.  RELEASE 00/47
C
C     AUTHORS ... GERALD T. ANTHONY, MAURICE G. COX
C                 J. GEOFFREY HAYES AND MICHAEL A. SINGER.
C     NATIONAL PHYSICAL LABORATORY, TEDDINGTON,
C     MIDDLESEX TW11 OLW, ENGLAND
C
C     *******************************************************
C
C     E01AEY.  A ROUTINE TO APPROXIMATE AND THEN REFINE AN
C     INTERPOLATING POLYNOMIAL  Q(X)  OR A ZEROIZING
C     POLYNOMIAL  Q0(X)  IN ITS CHEBYSHEV REPRESENTATION
C
C     INPUT PARAMETERS
C        WITHQ0   TRUE IF ZEROIZING POLYNOMIAL, ELSE FALSE
C        M        THE NUMBER OF DISTINCT DATA POINTS.
C        X        ARRAY CONTAINING THE DISTINCT X-VALUES,
C                    NORMALIZED IF NECESSARY TO (-1, 1).
C        XMIN,
C        XMAX     LOWER AND UPPER ENDPOINTS OF INTERVAL
C      * Y        ARRAY CONTAINING VALUES AND DERIVATIVES OF
C                    THE DEPENDENT VARIABLE.
C        IP       ARRAY SPECIFYING THE HIGHEST ORDER OF
C                    DERIVATIVE AT EACH X-VALUE.
C        IMAX     ONE MORE THAN THE LARGEST ELEMENT OF THE
C                    ARRAY IP.
C        N        NUMBER OF INTERPOLATING CONDITIONS.
C                    N = M + IP(1) + IP(2) + ... + IP(M).
C        NP1      VALUE OF  N + 1
C        ITMIN,
C        ITMAX    THE LOWER AND UPPER LIMITS ON THE ITERATIVE
C                    PROCESS.
C
C     OUTPUT (AND ASSOCIATED DIMENSION) PARAMETERS
C        A        CHEBYSHEV COEFFICIENTS OF POLYNOMIAL
C        LA       DIMENSION OF  A.
C                    .GE. N  IF INTERPOLATING POLYNOMIAL
C                    .GE. N + 1  IF ZEROIZING POLYNOMIAL
C        RES      RESIDUALS OF POLYNOMIAL
C        PMAX     LARGEST PERFORMANCE INDEX
C        PINDEX   PERFORMANCE INDICES
C        NIT      NUMBER OF ITERATIONS TAKEN
C
C     WORKSPACE PARAMETERS
C        ATRIAL   TRIAL VALUES OF THE CHEBYSHEV COEFFICIENTS
C        PTRIAL   PERFORMANCE INDICES CORRESPONDING TO  ATRIAL
C      * FTAU     SCALED VALUES OF  Y.  IF  Y(I)  IS THE
C                    VALUE OF AN  R-TH  DERIVATIVE, THEN
C                    ((XMAX - XMIN)/2)**R/(FACTORIAL R)
C                    TIMES  Y(I)  IS THE VALUE OF  FTAU(I)
C      * C        COEFFICIENTS IN NEWTON FORM OF POLYNOMIAL
C      * D        INTERMEDIATE DIVIDED DIFFERENCE VALUES
C      **W        VALUES OF CORRECTION POLYNOMIAL AT
C                    CHEBYSHEV EXTREMA
C        ADIF     CHEBYSHEV COEFFICIENTS OF A DERIVATIVE OF
C                    AN APPROXIMATION TO THE POLYNOMIAL
C        DA       CHEBYSHEV COEFFICIENTS OF A
C                    CORRECTION POLYNOMIAL
C        RNM      RESIDUAL NORMS CORRESPONDING TO  A
C        RTRLNM   RESIDUAL NORMS CORRESPONDING TO  ATRIAL
C      * LOCX     POINTERS TO X-VALUES IN CONSTRUCTING
C                    NEWTON FORM OF POLYNOMIAL
C      * LOCY     POINTERS TO Y-VALUES CORRESPONDING TO X-VALUES
C
C     FAILURE INDICATOR PARAMETER
C        IFAIL    FAILURE INDICATOR
C                    0 - SUCCESSFUL TERMINATION
C                    1 - ITERATION LIMIT EXCEEDED
C                    2 - ITERATION DIVERGENT
C
C           NOTES.  (1) THE ELEMENTS OF THE ARRAYS MARKED  *  ARE
C                       NOT ACCESSED IF  WITHQ0  IS  TRUE.
C                   (2) THE ELEMENTS OF THE ARRAY MARKED  **  IS
C                       NOT ACCESSED IF  WITHQ0  IS  FALSE.
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION  PMAX, XMAX, XMIN
      INTEGER           IFAIL, IMAX, ITMAX, ITMIN, LA, M, N, NIT, NP1
      LOGICAL           WITHQ0
C     .. Array Arguments ..
      DOUBLE PRECISION  A(LA), ADIF(LA), ATRIAL(LA), C(N), D(N), DA(LA),
     *                  FTAU(N), PINDEX(IMAX), PTRIAL(IMAX), RES(N),
     *                  RNM(IMAX), RTRLNM(IMAX), W(NP1), X(M), Y(N)
      INTEGER           IP(M), LOCX(M), LOCY(M)
C     .. Local Scalars ..
      DOUBLE PRECISION  AMAX, ATRNRM, DANRM, HALF, ONE, PMXTRL, SCALE,
     *                  SXTEEN, ZERO
      INTEGER           I, IERROR, IT, ITEMP, ITMXP1, ITP1, L, NFREF,
     *                  NPILT1, NTERMS
      LOGICAL           IMPROV, WITHPI, ZERODA
C     .. External Subroutines ..
      EXTERNAL          E01AEU, E01AEV, E01AEZ
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, LOG
C     .. Data statements ..
      DATA              ZERO, HALF, ONE, SXTEEN/0.0D+0, 0.5D+0, 1.0D+0,
     *                  16.0D+0/
C     .. Executable Statements ..
      IERROR = 0
C
C     NUMBER OF TERMS IN POLYNOMIAL
C
      NTERMS = N
      IF (WITHQ0) NTERMS = NTERMS + 1
C
C     NUMBER OF PERFORMANCE INDICES LESS THAN ONE
C
      NPILT1 = 0
C
C     INDICATE THAT PERFORMANCE INDICES ARE TO BE PRODUCED
C     ONLY IF AN IMPROVEMENT IS OBTAINED
C
      WITHPI = .FALSE.
C
C     INDICATE THAT THE FINE REFINEMENT STAGE HAS NOT YET STARTED
C
      NFREF = -2
C
C     SET RESIDUALS INITIALLY EQUAL TO SPECIFIED Y-VALUES
C     (IF  Q(X)  REQUIRED) OR ZERO (IF  Q0(X)  REQUIRED)
C
      DO 20 I = 1, N
         IF ( .NOT. WITHQ0) RES(I) = Y(I)
         IF (WITHQ0) RES(I) = ZERO
   20 CONTINUE
C
C     INITIALIZE TRIAL CHEBYSHEV COEFFICIENTS
C
      DO 40 I = 1, NTERMS
         ATRIAL(I) = ZERO
   40 CONTINUE
C
C     COMMENCE ITERATIVE REFINEMENT
C
      ITMXP1 = ITMAX + 1
      DO 360 ITP1 = 1, ITMXP1
C
C        IT  IS THE ACTUAL ITERATION NUMBER,  IT = 0
C        CORRESPONDING TO THE FIRST ESTIMATE OF THE POLYNOMIAL
C
         IT = ITP1 - 1
C
C        DETERMINE CHEBYSHEV COEFFICIENTS  DA  OF POLYNOMIAL
C        APPROXIMATELY SATISFYING THE CONDITIONS IN  RES
C
         IF (WITHQ0 .AND. IT.EQ.0) CALL E01AEU(M,X,IP,NP1,DA,W)
         IF ( .NOT. WITHQ0 .OR. (WITHQ0 .AND. IT.GT.0))
     *       CALL E01AEZ(M,XMIN,XMAX,X,RES,IP,N,DA,LOCX,LOCY,FTAU,D,C)
C
C        SKIP TEST FOR DIVERGENCE IF ON ZERO-TH ITERATION
C
         IF (IT.EQ.0) GO TO 100
C
C        DETERMINE THE NORMS OF  DA  AND (THE PREVIOUS)  ATRIAL
C
         DANRM = HALF*ABS(DA(1))
         ATRNRM = HALF*ABS(ATRIAL(1))
         IF (N.EQ.1) GO TO 80
         DO 60 I = 2, N
            DANRM = DANRM + ABS(DA(I))
            ATRNRM = ATRNRM + ABS(ATRIAL(I))
   60    CONTINUE
   80    IF (WITHQ0) ATRNRM = ATRNRM + ABS(ATRIAL(NP1))
C
C        ASSUME DIVERGENCE IF THE NORM OF  DA  IS NOT
C        LESS THAN THAT OF  ATRIAL  ...
C
         IF (DANRM.GE.ATRNRM) IERROR = 2
         IF (DANRM.GE.ATRNRM) GO TO 380
C
C        ... OTHERWISE DETERMINE NEW TRIAL APPROXIMATION
C
  100    ZERODA = .TRUE.
         DO 120 I = 1, N
            ATRIAL(I) = ATRIAL(I) + DA(I)
            IF (DA(I).NE.ZERO) ZERODA = .FALSE.
  120    CONTINUE
         IF ( .NOT. (WITHQ0 .AND. IT.EQ.0)) GO TO 140
         ATRIAL(NP1) = DA(NP1)
         IF (DA(NP1).NE.ZERO) ZERODA = .FALSE.
C
C        DETERMINE RESIDUALS, PERFORMANCE INDICES AND
C        LARGEST PERFORMANCE INDEX CORRESPONDING TO
C        TRIAL COEFFICIENTS  ATRIAL
C
  140    CALL E01AEV(WITHQ0,WITHPI,M,XMIN,XMAX,X,N,Y,IP,IMAX,ATRIAL,LA,
     *               IT,RNM,RTRLNM,IMPROV,ADIF,RES,PMXTRL,PTRIAL)
C
C        SET DUMMY, NON-ZERO, VALUE OF  PMXTRL  IF NO
C        IMPROVEMENT, OTHERWISE IT IS UNDEFINED
C
         IF ( .NOT. IMPROV) PMXTRL = SXTEEN
C
C        IF ON FIRST ITERATION, OR IF THE LARGEST PERFORMANCE
C        INDEX IS ZERO, OR IF ALL COMPONENTS OF  DA  ARE ZERO,
C        TAKE THE TRIAL SET OF COEFFICIENTS AND PERFORMANCE
C        INDICES AS THE BEST (SO FAR)
C
         IF ( .NOT. (IT.EQ.0 .OR. PMXTRL.EQ.ZERO .OR. ZERODA))
     *       GO TO 200
         DO 160 I = 1, NTERMS
            A(I) = ATRIAL(I)
  160    CONTINUE
         DO 180 L = 1, IMAX
            RNM(L) = RTRLNM(L)
            PINDEX(L) = PTRIAL(L)
  180    CONTINUE
         PMAX = PMXTRL
C
C        FINISH IF LARGEST PERFORMANCE INDEX IS ZERO
C        OR IF ALL COMPONENTS OF  DA  ARE ZERO
C        (I.E. NO FURTHER IMPROVEMENT IS POSSIBLE)
C
  200    IF (PMXTRL.EQ.ZERO .OR. ZERODA) GO TO 380
C
C        INDICATE WHETHER THE FINE REFINEMENT STAGE HAS COMMENCED
C        (I.E. FOR THE FIRST TIME ALL PERFORMANCE INDICES ARE
C        LESS THAN ONE)
C
         IF (NFREF.EQ.-2 .AND. PMXTRL.LT.ONE) NFREF = -1
C
C        BRANCH ACCORDING TO WHETHER THE PROCESS IS IN THE
C        FINE REFINEMENT STAGE  (NFREF .GE. 0)  OR NOT
C        (NFREF .EQ. -1)
C
         IF (NFREF.GE.-1) GO TO 280
C
C        THE PROCESS IS IN THE COURSE REFINEMENT PHASE.
C        UPDATE THE COEFFICIENTS AND THE CORRESPONDING
C        NORMS AND PERFORMANCE INDICES IF
C           (I)  THERE HAS BEEN AN IMPROVEMENT IN (AT
C                LEAST) ONE OF THE RESIDUAL NORMS, AND
C           (II) THE NUMBER OF PERFORMANCE INDICES
C                THAT ARE LESS THAN ONE HAS NOT
C                INCREASED COMPARED WITH THOSE OF
C                THE BEST POLYNOMIAL SO FAR.
C
         IF ( .NOT. IMPROV) GO TO 360
         ITEMP = 0
         DO 220 L = 1, IMAX
            IF (PTRIAL(L).LT.ONE) ITEMP = ITEMP + 1
  220    CONTINUE
         IF (ITEMP.LT.NPILT1) GO TO 360
         NPILT1 = ITEMP
         DO 240 I = 1, NTERMS
            A(I) = ATRIAL(I)
  240    CONTINUE
         DO 260 L = 1, IMAX
            RNM(L) = RTRLNM(L)
            PINDEX(L) = PTRIAL(L)
  260    CONTINUE
         PMAX = PMXTRL
         GO TO 360
C
C        THE PROCESS IS IN THE FINE REFINEMENT PHASE.
C        UPDATE THE COEFFICIENTS AND THE CORRESPONDING
C        NORMS AND PERFORMANCE INDICES IF
C           (I)  THERE HAS BEEN AN IMPROVEMENT IN (AT
C                LEAST) ONE OF THE RESIDUAL NORMS, AND
C           (II) THE LARGEST PERFORMANCE INDEX IS LESS
C                THAN THE LARGEST OF THAT OF THE BEST
C                POLYNOMIAL SO FAR.
C        INCREMENT THE NUMBER OF FINE REFINEMENTS (THE NUMBER
C        OF REFINEMENTS SINCE THE FIRST OCCASION WHEN ALL
C        PERFORMANCE INDICES WERE LESS THAN UNITY), EXITING
C        IF AS MANY AS  ITMIN  FINE REFINEMENTS HAVE BEEN
C        PERFORMED
C
  280    IF ( .NOT. IMPROV) GO TO 340
         IF (PMXTRL.GE.PMAX) GO TO 340
         DO 300 I = 1, NTERMS
            A(I) = ATRIAL(I)
  300    CONTINUE
         DO 320 L = 1, IMAX
            RNM(L) = RTRLNM(L)
            PINDEX(L) = PTRIAL(L)
  320    CONTINUE
         PMAX = PMXTRL
  340    NFREF = NFREF + 1
         IF (NFREF.GE.ITMIN) GO TO 380
  360 CONTINUE
C
C     THE PROCESS HAS NOT SUCCEEDED IN REDUCING ALL THE
C     PERFORMANCE INDICES TO LESS THAN UNITY
C
      IERROR = 1
      IT = ITMAX
C
C     NUMBER OF ITERATIONS ACTUALLY TAKEN
C
  380 NIT = IT
      IF ( .NOT. WITHQ0) GO TO 440
C
C     IN THE CASE OF  Q0(X),  SCALE ITS COEFFICIENTS BY
C     AN INTEGRAL POWER OF  16  SUCH THAT THE LARGEST
C     COEFFICIENT IS OF ORDER UNITY
C
      AMAX = ZERO
      DO 400 I = 1, NP1
         IF (ABS(A(I)).GT.AMAX) AMAX = ABS(A(I))
  400 CONTINUE
      IF (AMAX.EQ.ZERO) GO TO 440
      I = LOG(AMAX)/LOG(SXTEEN)
      SCALE = SXTEEN**(-I)
      DO 420 I = 1, NP1
         A(I) = SCALE*A(I)
  420 CONTINUE
C
C     RETURN RESIDUALS CORRESPONDING TO SELECTED COEFFICIENTS
C
  440 WITHPI = .TRUE.
      CALL E01AEV(WITHQ0,WITHPI,M,XMIN,XMAX,X,N,Y,IP,IMAX,A,LA,IT,RNM,
     *            RTRLNM,IMPROV,ADIF,RES,PMAX,PINDEX)
      IFAIL = IERROR
      RETURN
C
C     END E01AEY
C
      END

      SUBROUTINE E01AEZ(M,XMIN,XMAX,X,Y,IP,N,A,LOCX,LOCY,FTAU,D,C)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C
C     *******************************************************
C
C     NPL ALGORITHMS LIBRARY ROUTINE QPOLY
C
C     CREATED 02 05 80.  UPDATED 14 05 80.  RELEASE 00/09
C
C     AUTHORS ... GERALD T. ANTHONY, MAURICE G. COX
C                 J. GEOFFREY HAYES AND MICHAEL A. SINGER.
C     NATIONAL PHYSICAL LABORATORY, TEDDINGTON,
C     MIDDLESEX TW11 OLW, ENGLAND.
C
C     *******************************************************
C
C     E01AEZ. AN ALGORITHM TO DETERMINE THE CHEBYSHEV SERIES
C     REPRESENTATION OF A POLYNOMIAL INTERPOLANT  Q(X)  TO
C     ARBITRARY DATA POINTS WHERE DERIVATIVE INFORMATION MAY
C     BE GIVEN.
C
C     INPUT PARAMETERS
C        M        NUMBER OF DISTINCT X-VALUES.
C        XMIN,
C        XMAX     LOWER AND UPPER ENDPOINTS OF INTERVAL
C        X        INDEPENDENT VARIABLE VALUES,
C                    NORMALIZED TO  (-1, 1)
C        Y        VALUES AND DERIVATIVES OF DEPENDENT VARIABLE
C        IP       HIGHEST ORDER OF DERIVATIVE AT EACH X-VALUE
C        N        NUMBER OF INTERPOLATING CONDITIONS.
C                    N = M + IP(1) + IP(2) + ... + IP(M).
C
C     OUTPUT PARAMETERS
C        A        CHEBYSHEV COEFFICIENTS OF  Q(X)
C
C     WORKSPACE PARAMETERS
C        LOCX     POINTERS TO X-VALUES IN CONSTRUCTING
C                    NEWTON FORM OF POLYNOMIAL
C        LOCY     POINTERS TO Y-VALUES CORRESPONDING TO X-VALUES
C        FTAU     SCALED VALUES OF  Y.  IF  Y(I)  IS THE
C                    VALUE OF AN  R-TH  DERIVATIVE, THEN
C                    ((XMAX - XMIN)/2)**R/(FACTORIAL R)
C                    TIMES  Y(I)  IS THE VALUE OF  FTAU(I)
C        D        INTERMEDIATE DIVIDED DIFFERENCE VALUES
C        C        NEWTON COEFFICIENTS OF  Q(X)
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION  XMAX, XMIN
      INTEGER           M, N
C     .. Array Arguments ..
      DOUBLE PRECISION  A(N), C(N), D(N), FTAU(N), X(M), Y(N)
      INTEGER           IP(M), LOCX(M), LOCY(M)
C     .. Local Scalars ..
      DOUBLE PRECISION  CMIN, CNEW, FACTOR, ONE, PI, RI, RJ, S, SCALE,
     *                  SFAC, TWO, V, XCH, ZERO
      INTEGER           I, I2, IC, ICMIN, IFAIL, IFTAU, IP1, ISAVE, IY,
     *                  J, JMAX, K, KREV, L, LMAX, LOCXI, LOCXJ, LOCXK,
     *                  LOCYI, NC
C     .. External Functions ..
      DOUBLE PRECISION  X01AAF
      EXTERNAL          X01AAF
C     .. External Subroutines ..
      EXTERNAL          E01AEX, E02AFF
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, SIN
C     .. Data statements ..
      DATA              ZERO, ONE, TWO/0.0D+0, 1.0D+0, 2.0D+0/
C     .. Executable Statements ..
      PI = X01AAF(ZERO)
      SCALE = (XMAX-XMIN)/TWO
C
C     INITIALIZE X- AND Y-POINTERS
C
      IY = 0
      DO 40 I = 1, M
         LOCX(I) = I
         IY = IY + 1
         LOCY(I) = IY
         FTAU(IY) = Y(IY)
         JMAX = IP(I)
         IF (JMAX.EQ.0) GO TO 40
C
C        FORM THE SCALED DERIVATIVES, I.E. AN  R-TH  DERIVATIVE
C        VALUE IS DIVIDED BY  FACTORIAL R  AND MULTIPLIED
C        BY THE  R-TH  POWER OF  (XMAX - XMIN)/2
C
         SFAC = ONE
         DO 20 J = 1, JMAX
            IY = IY + 1
            RJ = J
            SFAC = SFAC*SCALE/RJ
            FTAU(IY) = Y(IY)*SFAC
   20    CONTINUE
   40 CONTINUE
C
C     FORM SUCCESSIVE UPWARD SLOPING DIAGONALS OF
C     THE DIVIDED DIFFERENCE TABLE
C
      NC = 0
      DO 120 J = 1, M
C
C        CHOOSE EACH X-VALUE IN TURN TO MAKE THE CORRESPONDING
C        NEWTON COEFFICIENT AS SMALL IN MAGNITUDE AS POSSIBLE
C
         DO 80 I = J, M
            LOCXI = LOCX(I)
            LOCYI = LOCY(LOCXI)
            CALL E01AEX(M,X,IP,N,LOCX,C,NC,X(LOCXI),J,FTAU(LOCYI)
     *                  ,1,CNEW,D)
            IF (I.EQ.J) GO TO 60
            IF (ABS(CNEW).GE.ABS(CMIN)) GO TO 80
   60       CMIN = CNEW
            ICMIN = I
   80    CONTINUE
         C(NC+1) = CMIN
         ISAVE = LOCX(J)
         LOCXJ = LOCX(ICMIN)
         LOCX(ICMIN) = ISAVE
         LOCX(J) = LOCXJ
C
C        CALCULATE THE RESULTING NEWTON COEFFICIENT (I.E.
C        REPEAT THE ABOVE COMPUTATION, BUT ONLY IN THE CASE
C        LEADING TO THE SMALLEST NEW COEFFICIENT)
C
         IFTAU = LOCY(LOCXJ) - 1
         IP1 = IP(LOCXJ) + 1
         DO 100 I = 1, IP1
            IFTAU = IFTAU + 1
            CALL E01AEX(M,X,IP,N,LOCX,C,NC,X(LOCXJ),J,FTAU(IFTAU)
     *                  ,I,C(NC+1),D)
            NC = NC + 1
            IF (NC.EQ.N) GO TO 140
  100    CONTINUE
  120 CONTINUE
C
C     EVALUATE  Q(X)  (FROM ITS NEWTON FORM) AT THE EXTREMA
C     OF THE CHEBYSHEV POLYNOMIAL OF DEGREE  N - 1  ...
C
  140 FACTOR = 2*N - 2
      FACTOR = PI/FACTOR
      I2 = N + 1
      DO 200 I = 1, N
         I2 = I2 - 2
         RI = I2
         XCH = SIN(FACTOR*RI)
         S = C(N)
         IC = N
         K = M + 1
         DO 180 KREV = 1, M
            K = K - 1
            LOCXK = LOCX(K)
            LMAX = IP(LOCXK) + 1
            IF (K.EQ.M) LMAX = LMAX - 1
            IF (LMAX.LE.0) GO TO 180
            V = XCH - X(LOCXK)
            DO 160 L = 1, LMAX
               IC = IC - 1
               S = S*V + C(IC)
  160       CONTINUE
  180    CONTINUE
         D(I) = S
  200 CONTINUE
C
C     ... IN ORDER TO DETERMINE THE COEFFICIENTS IN ITS
C     CHEBYSHEV REPRESENTATION
C
      CALL E02AFF(N,D,A,IFAIL)
      RETURN
C
C     END E01AEZ
C
      END
